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PROJECT SUMMARY 


The purpose of this study is to establish the technical ground for modeling 
the physics of laser powered pulse detonation phenomenon. The principle of the 
laser power propulsion is that when high-powered laser is focused at a small 
area near the surface of a thruster, the intense energy causes the electrical 
breakdown of the working fluid (e.g. air) and forming high speed plasma (known 
as the inverse Bremsstrahlung, IB, effect). The intense heat and high pressure 
created in the plasma consequently causes the surrounding to heat up and 
expand until the thrust producing shock waves are formed. This complex process 
of gas ionization, increase in radiation absorption and the forming of plasma and 
shock waves are investigated in the development of the present numerical 
model. In the first phase of this study, laser light focusing, radiation absorption 
and shock wave propagation during the laser pulsed cycle are modeled. The 
model geometry and test conditions of known benchmark experiments such as 
those in Myrabo’s experiment are employed in the numerical model benchmark 
validation studies. The calculated performance data (e.g. coupling coefficients) 
are compared to the measured data. The final goal of this project is that a design 
tool will be available for the analysis and optimization of full-scale laser propelled 
flight vehicles. 
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1. METHOD OF APPROACH - A BUILDING BLOCK APPROACH 


An efficient design of laser propelled vehicles requires the detailed 
understanding of how the detonation waves behave with different thrust cavity 
geometry and flight conditions and the performance of the air intake 
configurations and arrangements. CFD methods with high-temperature 
thermodynamics and transient detonation wave capturing capabilities are very 
useful in providing information for the optimization of the laser lightcraft 
configurations. Compared to the other propulsion systems, modeling of laser 
propelled vehicles is very difficult because it involves many complicated physical 
phenomena such as optical breakdown, laser absorption, thermochemical and 
raidiative nonequilibrium, plasma dynamics, etc. The general approach of the 
present numerical model can be summarized in the following items: 

• Computational fluid dynamics method with time accurate integration schemes 
and high-temperature gas thermodynamics properties 

• Unstructured-mesh flow code for detonation wave prediction 

• Thermal non-equilibrium energy equations monitoring the translational and 
vibrational states 

• Modeling of the electron temperature transport and energy transfer 

• Plasma ignition model to start the laser absorption 

• Air breakdown/ionization chemistry model 

• Specular ray tracing method for laser light reflection and focusing 

• Modified discrete transfer radiation model for laser absorption and the inverse 
Bremsstrahlung effect 

• Laser light energy absorption model 

• Beam propagation and gasdynamics coupling 

• Effects of air plasma non-equilibrium thermal radiation model 

• Thrust integration in time for coupling coefficient calculations 
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2. COMPUTATIONAL FLUID DYNAMICS MODEL 


The current computational model is developed based on an unstructured 
computational fluid dynamics (CFD) model, the UNIC-UNS code developed at 
Engineering Sciences, Inc., and a ray tracing radiation model for laser light 
reflection, focusing and absorption modeling. The underlying CFD flow solver is a 
general unified (all speed) solution method employed to solve the Navier-Stokes 
equations (Ref. 1). New models involved in this project include: 

(1) high-temperature thermal dynamics properties based on the latest 
database released in Ref. 2; 

(2) an electron temperature equation with elastic collision energy transfer 
described in Ref. 3; 

(3) a finite-rate air chemistry model of Park (Ref. 4); 

(4) a laser ray tracing model with discrete transfer method in solving the 
radiative transfer equation; 

(5) the plasma radiation model using the LORAN code (Ref. 5); and 

(6) a spark ignition model to start the laser light absorption process and 
chemical reactions. 


2.1 Governing Equations 

The continuity, Navier-Stokes, energy (total enthalpy) and electron 
temperature equations, can be written in a Cartesian tensor form: 
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where p is the fluid density, Uj is the t h Cartesian component of the velocity, p is 
the static pressure, p is the fluid viscosity, P r is the Prandtl number, H is the gas 
total enthalpy and V stands for the sum of velocity squared. In Eq. (4), k b , n e , T e , 
Xe, Q r and Qec are the Boltzmann’s constant, electron number density, electron 
temperature, electron thermal conductivity, radiative heat source from laser 
absorption and radiative transfer, and the energy transfer due to eletron/particle 
elastic collisions, respectively. The shear stress fjcan be expressed as: 


=M 


dx- 3 dx k * j 




The species conservation equation is expressed as: 
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where Yj is the f h species mass fraction, D is the mass diffusivity, which can be 
written as viscosity divided by the Schmidt number, o Y , and d) t is the chemical 
reaction rate for species / respectively. 


2.2 Thermal Non-Equilibrium Energy Equations 

For high temperature gas flows, thermal non-equilibrium state may be 
important. In Landau and Teller’s derivation, a master equation is employed to 
describe the evolution of the population of quantum level N|. This master 
equation is written as: 

dN Lw 

K^jNr. i = 0,1,2, 

ai 0 ;-0 
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Results from the quantum mechanical solution of the harmonic oscillator are 
used to relate the various quantum transition rates to one another, and then the 
master equation may be summed over all quantum states to arrive at the Landau- 
Teller equation: 

£pe ! __j_{ £ZLL 

D, 3x\ K '9x i )* P t LT 


where p, e v , e v eq and %n represent the gas density, vibrational energy, effective 
(equilibrium) vibrational energy and the vibrational-translational relaxation time 
scale respectively. An empirical expression (to be discussed in the next section) 
is used to model the Landau-Teller relaxation time scale. 

To solve this vibrational energy equation, a new subroutine is created to 
calculate the matrix coefficients (the left-hand-side of the equation) of the 
transport equation and the source term (the right-hand-side of the equation). The 
source term is further linearized to result in an explicit term and an implicit term. 
This treatment is important for an unconditionally stable solution of the equation. 
That is, 

e eq (T) 

Explicit source term = p-- ~ — 

*LT 

O 

Implicit source term = e v 

^LT 

The vibrational to translational energy transfer is computed by adding the 
diffusion term of the vibrational energy equation, 
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to the right-hand-side of the translational energy equation as a source term. 


2.3 Numerical Scheme 

The cell-centered scheme is employed that the volume surfaces are 
represented by the grid cell surfaces. The transport equations can be written in 
integral form as: 

•—j pfidtl + i F hdT =\ S 6 d£l 

where Q is the domain of interest and r denotes the surrounding surfaces; n is a 
unit normal vector of r in the outward direction. The flux function F contains the 
inviscid and the viscous flux vectors, 

F = pV(})-fyV<p 



Figure 1 . Unstructured control volume. 


For the face e between control volumes P and E, the diffusive flux can be 
approximated as: 


And, 


•n) t - 


. n- 
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where V<p e is interpolated from the neighbor cells E and P. 
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The finite volume formulation of flux integral can be evaluated by the 
summation of the flux vectors over each face, 

Lt-ar* X F u Ar , 

;=*</) 

where k(i) is a list of faces of cell i, Fy represents convection and diffusion fluxes 
through the interface between cell i and j, and AT, is the cell-face area. 

The convective flux is evaluated through the upwind-cell quantity by a linear 
reconstruction procedure to achieve second order accuracy: 




where the subscript u represents the upwind cell and y/ e is a limiter used to 
ensure that the reconstruction does not introduce local extrema. The limiter 
proposed by Barth is used here. Defining <p max = max{f u ,jj), <l> min =min{<i) u ,<p .) (and 


assuming <p° e is computed with ^e=1) the scalar y e associated with the gradient 
at cell u due to edge e is: 
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2.4 Solution Procedures 

A general implicit discretized time-marching scheme for the transport 
equations can be written as below, 
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where NB means the neighbor cells of cell P. The high order differencing term 
and cross diffusion term are treated using known quantities and retained in the 
source term and updated explicitly. 

A predictor/corrector solution algorithm is employed to provide coupling of the 
governing equations. 

The discretized finite-volume equations form a set of linear algebra equations, 
which are non-symmetric matrix system with arbitrary sparsity patterns. The 
preconditioned Bi-CGSTAB and GMRES(/7i) matrix solvers are used to efficiently 
solve the linear algebra equations. 


2.5 Finite-Rate Chemistry Model 

For gas-phase chemical reaction modeling, a general system of chemical 
reactions can be written in terms of its stoichiometric coefficients (vy and vy’ ) and 
the t h chemical species name (A//,) of the/* reaction as 

X v u M i = X v ‘> M < 

i i 

The net rate of change in the molar concentration of species / due to 
reactions j, X,y, and the species production rate can be written as: 
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The forward (Arrhenius law) and backward reaction rates for each reaction is 
given by: 

K fj =AjT B ’Exp {-jL) 

K bi = «L 
' Ke, 


where Key is the equilibrium coefficient 


Ke ) = (RT) Exp 


1 (V.. - V..)g. 
j_ E j_ u 

RT 


A point-implicit (operator splitting) method is employed to solve the chemistry 
system. 


2.6 Air Chemistry Model 

The air chemical kinetics model of Park (Ref. 4) is employed in the present 
study. This model is listed in Table 1. The light emission effects due to the 
radiative recombination reactions, in the last two reactions, are not modeled in 
the current study. This effect will be investigated in future study. Other chemistry 
models such as the one described in Ref. 6 were also tested. That model, 
however, caused very high stiffness in the high electron temperature regime. It is 
therefore not selected for the final model. The current chemistry model involves 
1 1 species and 51 reactions. The present point implicit finite-rate method have is 
very robust in providing solutions of this system. 
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Table 1. Air Chemistry Model 

REACT RATE: K = A*T**(+N)*EXP(-E/RT) 


I A N E/R Temperature (T) (EQUATIONS) 


1 

7 . 0000e21 

-1.60 

113200.0 

Tc 

N2 

+ 

N2 


= 

2 

N 

+ 

N2 

2 

7 . 0000e21 

-1.60 

113200.0 

Tc 

N2 

+ 

02 


= 

2 

N 

+ 

02 

3 

7 . 0000e21 

-1.60 

113200.0 

Tc 

N2 

+ 

NO 


= 

2 

N 

+ 

NO 

4 

7 . 0000e21 

-1.60 

113200.0 

Tc 

N2 

+ 

NO+ 


= 

2 

N 

+ 

NO+ 

5 

3 . 0000e22 

-1.60 

113200.0 

Tc 

N2 

+ 

N 


= 

2 

N 

+ 

N 

6 

3 . 0000e22 

-1.60 

113200.0 

Tc 

N2 

+ 

N+ 


= 

2 

N 

+ 

N+ 

7 

3 . 0000e22 

-1.60 

113200.0 

Tc 

N2 

+ 

O 


= 

2 

N 

+ 

0 

8 

3 . 0000e22 

-1.60 

113200.0 

Tc 

N2 

+ 

0+ 


= 

2 

N 

+ 

0+ 

9 

7 . 0000e21 

-1.60 

113200.0 

Tc 

N2 

+ 

N2 + 


= 

2 

N 

+ 

N2 + 

10 

7 . 0000e21 

-1.60 

113200.0 

Tc 

N2 

+ 

02 + 


= 

2 

N 

+ 

02 + 

11 

1 . 2000e25 

-1.60 

113200.0 

Tc 

N2 

+ 

E 



2 

N 

+ 

E 

12 

2 . 0000e21 

-1.50 

59500.0 

Tc 

02 

+ 

N2 



2 

0 

+ 

N2 

13 

2 . 0000e21 

-1.50 

59500.0 

Tc 

02 

+ 

02 


= 

2 

0 

+ 

02 

14 

2 . 0000e21 

-1.50 

59500.0 

Tc 

02 

+ 

NO 


= 

2 

0 

+ 

NO 

15 

2 . 0000e21 

-1.50 

59500.0 

Tc 

02 

+ 

NO+ 


= 

2 

0+ N0+ 

16 

1 . 0000e22 

-1.50 

59500.0 

Tc 

02 

+ 

N 



2 

0 

+ 

N 

17 

1 .0000e22 

-1.50 

59500.0 

Tc 

02 

+ 

N+ 


= 

2 

0 

+ 

N+ 

18 

1 .0000e22 

-1.50 

59500.0 

Tc 

02 

+ 

O 


= 

2 

0 

+ 

0 

19 

1 .0000e22 

-1.50 

59500.0 

TC 

02 

+ 

0+ 


= 

2 

0 

+ 

0+ 

20 

2 . 0000e21 

-1.50 

59500.0 

TC 

02 

+ 

N2 + 


= 

2 

0 

+ 

N2 + 

21 

2 . 0000e21 

-1.50 

59500.0 

TC 

02 

+ 

02 + 


= 

2 

0 

+ 

02 + 

22 

5 . OOOOelS 

0.00 

75500.0 

Tc 

NO 

+ 

N2 


= 

N 

+ 

0 

+ N2 

23 

5 . OOOOelS 

0.00 

75500.0 

Tc 

NO 

+ 

02 


= 

N 

+ 

0 

+ 02 

24 

1.1000el7 

0.00 

75500.0 

Tc 

NO 

+ 

NO 



N 

+ 

0 

+ NO 

25 

5 . OOOOelS 

0.00 

75500.0 

Tc 

NO 

+ 

NO+ 


= 

N 

+ 

0 

+ N0+ 

26 

1 .1000el7 

0.00 

75500.0 

Tc 

NO 

+ 

N 


= 

N 

+ 

0 

+ N 

27 

1 . 1000el7 

0.00 

75500.0 

Tc 

NO 

+ 

N+ 



N 

+ 

0+ N+ 

28 

1 .1000el7 

0.00 

75500.0 

Tc 

NO 

+ 

O 


= 

N 

+ 

0 

+ 0 

29 

1 . 1000el7 

0.00 

75500.0 

Tc 

NO 

+ 

0+ 


= 

N 

+ 

0 

+ 0+ 

30 

5 .OOOOelS 

0.00 

75500.0 

Tc 

NO 

+ 

N2 + 


= 

N 

+ 

0 

+ N2 + 

31 

5. OOOOelS 

0.00 

75500.0 

Tc 

NO 

+ 

02 + 


= 

N 

+ 

0 

+ 02 + 

32 

8 . 4000el2 

0.00 

19400.0 

Tg 

NO 

+ 

O 


= 

02 4 

■ N 

33 

6 . 4000el7 

-1.00 

38400.0 

Tg 

N2 

+ 

O 


= 

NO ^ 

■ N 

34 

8 . 8000e08 

1.00 

31900.0 

Tg 

N + 

( 



= 

NO+ 

+ 

E 

35 

7 . 1000e02 

2.70 

80600.0 

Tg 

2 O 



= 

02 + 

+ 

E 

36 

4 . 4000e07 

1.50 

67500.0 

Tg 

2 N 



= 

N2 + 

+ 

E 

37 

1 .0000el2 

0.50 

77200.0 

Tg 

NO+ - 

»- 0 


= 

N+ 4 

■ 02 

38 

1 .0000el2 

0.50 

12200.0 

Tg 

N+ 

+ 

N2 


= 

N2 + 

+ 

N 

39 

8 . 7000el3 

0.14 

28600.0 

Tg 

02+ - 

»- N 


= 

N+ 4 

■ 02 

40 

1 .4000e05 

1.90 

26600.0 

Tg 

NO 

+ 

0+ 


= 

N+ 4 

■ 02 

41 

9 . 9000el2 

0.00 

40700.0 

Tg 

02+ - 

f N2 


= 

N2 + 

+ 

02 

42 

4 . 0000el2 

-0.09 

18000.0 

Tg 

02+ * 

► O 



0+ 4 

■ 02 

43 

3 . 4000el3 

-1.08 

12800.0 

Tg 

NO+ - 

► N 



N2 * 

■ 0+ 

44 

2 . 4000el3 

0.41 

32600.0 

Tg 

NO+ * 

t- 02 


= 

02 + 

+ 

NO 

45 

7 . 2 000el2 

0.29 

48600.0 

Tg 

NO+ - 

► 0 


= 

02+ 

+ 

N 

46 

9 . lOOOell 

0.36 

22800.0 

Tg 

0+ 

+ 

N2 


= 

N2 + 

+ 

0 

47 

7 . 2 000el3 

0.00 

35500.0 

Tg 

NO+ - 

H N 


= 

0 

+ 

N2 + 

48 

3 . 9000e33 

-3.78 

158500.0 

Te 

O + 

E 


= 

0+ 4 

■ ] 

E + E 

49 

2 . 5000e34 

-3,82 

168600.0 

Te 

N + 

1 

7 


= 

N+ 4 

• ] 

E + E 

50 

1 . 0700ell 

-0.52 

0.0 

Te 

0+ 

+ 

E + 

E 

= 

0 

+ 

E 


51 

1 . 5200ell 

-0.48 

0.0 

Te 

N+ 

+ 

E + 

E 

= 

N 

+ 

E 



where Tg stands for the gas temperature, Te is the electron temperature and Tc denotes the 


geometric averaged temperature, i.e. t c = ^ T g T v • TV is the vibrational temperature. 
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3. RADIATIVE HEAT TRANSFER MODEL 


3.1 Advantages and Disadvantages of Discrete Transfer Method 

The discrete ordinate method (DOM) and discrete transfer method (DTM) are 
currently two of the most widely used methods for modeling radiation in 
participating media, and they have many similarities. In both of these two 
methods, the radiative transfer equation (RTE) is solved only along discrete 
directions that approximate angular intensity distributions. As a result, they are 
mathematically very simple and relatively easy to code; they are able to treat 
multi-dimensional problems with complicated geometries (including shadowing 
effects) with ease; they can account for spectral effects from gases or surfaces; 
they are compatible with numerical algorithms for solving other transport 
equations. The difference between the DOM and DTM lies in how the RTE is 
solved. The former uses the finite difference method to solve the RTE while the 
latter uses the ray tracing method to solve the RTE. 

During the pulsed laser plasma ignition process the major quantity of interest 
from laser radiation is the net radiative heat transfer along the laser beam, q r , 
and it enters the governing equations through the energy source term. The 
wavelength dependent expression for q r is given by 



4=1 


where the subscript k refers to a particular radiation band and N b designates the 
total number of bands used in describing the wavelength dependence of the 
absorptivity a and emissivity e; G denotes the incident radiation on the surface; 
E b represents the radiative emissive power defined by the Planck function. In the 
DOM and DTM, the quantity G can be further expressed as 


li 


( 6 ) 


*=1 


N b N b m=\,M Nj, m=\,M 

G = %G k =Z £c te =£ 1 /^ 

k=\ n n m < 0 


k = 1 <0 


where n represents the unit normal vector on the wall; £2 m refers to the m-th 
specific discrete direction; Atom is the corresponding solid angle for the m-th 
discrete direction; Gi™, Ikm represent the incident radiation and incident intensity 
along the m-th direction at the k-th band, respectively. The DOM is not suitable 
for modeling surface radiation because it has to model radiation transport 
through each gas control volume even though the gas is not radiatively 
participating, and it is subsequently expensive. The DTM has been also rarely 
used to model surface radiation solely due to its fundamental shortcoming, the 
ray effects. 

To illustrate the ray effects in DTM, consider a radiation incident on a surface 
cell on the bottom wall over a control angle for a 3D box as shown in Fig. 6. The 
center point of the incident surface cell is O. The control angle is intercepted by 
the top wall and the intercepted area is represented by the shadowed area 
ABCD. Each surface of the box is divided into many small surface cells, the area 
ABCD thus consists of several surface cells. Some cells are entirely located 
inside the ABCD while others are partially located in the ABCD. The temperature 
and radiative properties on each cell may be different from others. In the 
traditional approach of the DTM, radiation contribution from the area ABCD to the 
incident surface cell is assumed to come from only a single surface cell, for 
example, where the point P is located in Fig. 2. Unless a large amount of discrete 
directions are chosen, this treatment neglects the contribution from other surface 
cells, thus leads to a poor accuracy in the results. This is so called the ray 
effects. The ray effects also exist in the gas radiation modeling but they become 
less prominent due to the mitigation of gas radiation. 
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3.2 Modified Discrete Transfer Method 

In order to eliminate the ray effects, the traditional DTM is modified by 
considering radiation contribution from all surface cells in the ABCD. Thus Eq. (6) 
becomes, 

Nj, m=\,M N 

G = £ I £ 


where j refers to the j-th surface cell in the ABCD; l km j is the radiative intensity 
from the j-th cell along the m-th direction at the k-th band, and Aa> mj is the solid 
angle with which the point O is seen by the j-th surface cell. The modified DTM 
(MDTM) keeps all the advantages of DTM such as easy consideration of 
spectral, specular, shadowing effects, easy coding and easy coupling with other 
transport equations while it provides a high accuracy in the results. Calculation of 
Atom] requires an evaluation of all surface cell areas intercepted by a control 
angle on the wall, and subsequently the ray tracing time is increased compared 
to the DTM. However, this increase is of minor effects since the surface cell 
number only takes a small percentage of the total cell number and Acomj is 
evaluated only once in the transient simulation of practical systems. Extension of 
the MDTM to the gas radiation problems requires a significant increase of ray 
tracing calculation, thus is not encouraged. 

The dependent variable in the MDTM is the radiative intensity as seen in Eqs. 
(6) and (7) and its value is initially not known. Thus, an iterative solution 
procedure is required in the above modeling of surface radiation. That is, an 
initial intensity solution is assumed and it is used to calculate the incident 
radiation G, then the incident radiation G and the Planck function are applied on 
the radiative boundary condition to obtain the new intensity solution. With new 
intensity solution, the new incident radiation can be calculated again. This 
iterative procedure continues until the maximum difference in the incident 
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radiation between two consecutive iterations is less than a designated tolerance. 
For a diffuse wall, the radiative boundary condition is given as 


I km ~ £ k 



+ 


g-g*) 

n 


1 G„ 


«n„.<o 


( 8 ) 


For a specular wall, the radiative boundary condition becomes 


/ 


km 



(1 z£j> 

7t 


9 km ’ 


(9) 


In Eq. (8), the all arriving intensities have a contribution to the calculation of Ikm 
while in Eq. (9), only one arriving intensity, whose direction is aligned with the 
mirror direction of the m-th direction, has a contribution to the calculation of 1^. 


3.3 Surface Cell Areas Intercepted by a Control Angle 

As indicated previously, the major calculation in the implementation of the 
MDTM comes from the evaluation of the solid angles with which surface cells 
partially or entirely intercepted by a control angle are seen from a certain point. 
Considering that a solid angle is defined as the projection of the surface onto a 
plane normal to the direction vector, divided by the distance squared, the major 
calculation in the MDTM thus becomes the calculation of surface cell areas. For 
some surface cells, entire area must be calculated while for others only portion of 
area need to be calculated as seen in Fig. 2. The following paragraphs will 
discuss how the surface cell areas are evaluated in this study. 

The first step is to select an appropriate discretization strategy for the angular 
domain. Due to its easy division and wide applications, the azimuthal 
discretization strategy as shown in Fig. 3 is chosen in this study. In this strategy, 
the 4k angular domain at a volume cell is divided into a finite number of discrete, 
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nonoverlapping solid angles along the azimuthal and polar directions. The angle 
increments along these two directions are usually uniform. For a surface 
radiation problem, only surface cells are taken into account and their angular 
domain is half of that for a volume cell. In Fig. 3, the shadowed area represents 
the projection of a control angle on a hemisphere of unit radius and it contains 
four vertex points. In this study, the direction vectors that connect the point O on 
the surface and four vertex points are called the vertex vector. A ray travelling 
along a vertex vector is called the vertex ray. Four vertex rays uniquely decide a 
control angle. 

The second step is to find the intersection points between the four vertex rays 
from a control angle and a wall as shown as the points A, B, C and D in Fig. 2. 
For the simulation in manufacturing and materials processing, the gas region 
usually needs to be discretized due to the involvement of other transport 
phenomena. With discretization of the gas region, calculation of the intersection 
point between a vertex ray and a wall becomes easy and straightforward. A 
vertex ray emitted from a surface cell is just traced through the adjacent volume 
cells until it is intercepted by a wall. Such a ray tracing calculation is very efficient 
and accurate. 

The third step is to determine the penetrating points of all surface cell edges 
passing through the four planes as shown in Fig. 2 as OAB, OAD, ODC, and 
OBC that form a control angle. These planes are defined by the two adjacent 
vertex vectors and limited between the corresponding vertex rays. Calculation of 
this step starts from one of the intersection points between the vertex rays and a 
wall, for example, the point A. The surface cell index number where the point A is 
located must be found at first. Then, each edge of this cell is examined to see if it 
penetrates any of four planes forming the given control angle. If one of these 
edges is found, the penetrating point must be determined. After that, the search 
goes to the adjacent surface cell with the common edge. This procedure 
continues until all the penetrating points are found. 
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The fourth step is to determine the vertexes of surface cells which are located 
inside the specified control angle. This step starts from a surface cell which has 
at least an edge penetrating one of the four planes forming the control angel. 
Each vertex of this cell is examined to see if it is located inside the control angle. 
If one of these vertexes is found, it is then saved and the search moves to the 
adjacent cells with the common vertex. This procedure repeats until a surface 
cell is found which has no vertex located inside the control angle. 

The fifth step is to calculate the areas of the surface cells contained inside the 
control angle. In this step, all penetrating points determined in Step (3) and 
vertexes determined in Step (4) are sorted and regrouped. Those points which 
are attached to the same surface cell are put together, and the area of the 
polygon formed by these points can then determined easily. 

In the above five steps, the Step (3) is the most computationally extensive. 
This is especially true for a 2D axisymmetric geometry. In such a case, the edges 
of a surface cell may be no longer straight lines as in 2D planar and 3D 
geometries. Thus, much more calculation has to be involved to find a valid 
penetrating point. Development of an accurate and efficient method for 
calculating the surface cell areas intercepted by the control angles represents 
one of the most important tasks in the implementation of the MDTM. The above 
method for calculating the surface cell areas has a room for improvement by 
using some advanced ray tracing techniques developed in computer graphics. 
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Figure 2. Concept for the modified discrete transfer method 



Figure 3. Azimuthal discretization strategy 


4. Non-equilibrium Radiative Heat Transfer 


The non-equilibrium radiation model needed for high temperature plasma is 
implemented based on the LORAN code of NASA Langley Research Center. The 
code is acquired from Langley, several published cases were tested before the 
implementation is started. After some adjustment of the input parameters, the 
code reproduced the results in the literature for high-temperature gas radiation. 

Implementation and coupling of the LORAN code with UNIC-UNS/ GRASP- 
UNS need further investigation due to the large amount of computational effort 
would be needed if it is solved for every grid point. The LORAN code solves 
thousands of bands in the radiative transfer model. More efficient method would 
be required to make this model computationally feasible. 


The general consideration of the non-equilibrium radiation effects and model 
adapted are summarized in the following. Conclusion of a preliminary study 
about the significance of its effects is also given below. 


1. Neglecting transients and assuming a non-scattering medium, the complete 
radiative transfer equation (RTE) becomes 


dlyis.a) 

ds 


+tcJJs,Q) = j' a (s) 


where / w represents spectral radiation intensity, and and f m denote the 
absorption coefficient and the emission coefficient, respectively. 


2. In equilibrium gases, the electronic energy-level populations are determined 
as a function of a uniquely defined equilibrium temperature according to a 

Boltzmann distribution, and K m and j e m are related according to Kirchhoff’s 
law as tiis) = Kj b js ) , where /*($) is blackbody intensity determined by 
Planck function. 
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In non-equilibrium flow, however, all these simple relations no longer apply 
and the non-equilibrium absorption and emission coefficients must be 
determined through the non-equilibrium populations of each energy level and 
transitions of various energy levels. 

3. Plasma air includes the following atomic and molecular species: O, N, NO, 
N + , 0 + , N 2 , 0 2 , etc. which all contribute to the non-equilibrium radiative heat 

transfer. To determine K m and j e a for these species, the following four 
radiative transitions must be taken into account: 

> Atomic Line Transitions 

> Atomic Bound-free Transitions 

> Atomic Free-free Transitions 

> Molecular Transitions 

4. Currently, there are two codes available, which provide detailed information 
on these transitions. One is NEQAIR developed by Park and another is 
LORAN developed by Hartung. Because these codes involve large database, 
they are rarely applied in practical problems. To overcome disadvantage of 
these code, some simple models have been developed such as Park’s PRG 
model etc. 

5. With the determination of K a and £, the RTE can be solved by either 

deterministic or stochastic approach. Unlike the most CFD equations, RTE is 
integral differential equation and numerical treatment is, thus, different from 
CFD approach. Currently, there are several methods available for solving 
non-equilibrium radiative heat transfer, which include Monte Carlo method, P- 
1 method, S-N method, quadromoment method, etc. The Monte Carlo is 
accurate but too costly and time-consuming for practical applications. The P-1 
and quadromoment methods are not applicable for optically thin medium. The 
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S-N method can provide accurate results for all optical ranges if the discrete 
direction number is big enough. 

6. In the present project, the S-N method and the finite volume method were 
used in short duration test runs. It is found that the non-equilibrium radiation 
model is very computational intensive. It is typically more than ten times 
slower than the model without non-equilibrium radiation model. Also, 
preliminary studies have shown that its contribution to the overall energy 
reduction (cooling effect) is only less than 10 percent. Therefore, this model is 
not used in the final performance prediction runs. 

4.1 Non-equilibrium Air Plasma Radiation Analysis 

The current modeling efforts for radiative heat transfer in laser induced 
plasmas involve many assumptions. Some of these assumptions are usually not 
valid and thus cause considerably errors in the modeling results. In the present 
study, efforts have been made to employ less assumptions and more accurate 
methods and models. Treatment of plasma radiative heat transfer involves 
solving the radiative transfer equation (RTE). Neglecting transients and assuming 
a noscattering medium, the complete RTE becomes 


ds 


+ *Ja>(s, &) = &(*) 


where /^represents spectral radiation intensity, and K m and y* denote the 
plasma spectral absorption coefficient and the emission coefficient, respectively. 

In equilibrium gases, the electronic energy-level populations are determined 
as a function of a uniquely defined equilibrium temperature according to a 
Boltzmann distribution, and K m and y' are related according to Kirchhoff’s law as 

j e a) (s) = /cjl(s ) , where l b a (s ) is blackbody intensity determined by Planck 
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function. In non-equilibrium flow, however, all these simple relations no longer 
apply and the non-equilibrium absorption and emission coefficients must be 
determined by the non-equilibrium populations of each energy level and 
transitions of various energy levels. 

Air plasma includes following atomic and molecular species: O, N, NO, N + , 
0\ NO + , N 2 , 0 2 , etc. which all contribute to the non-equilibrium radiative heat 
transfer. To determine K m and j e w for these species, the following four radiative 

transitions must be taken into account: atomic line transitions, atomic bound-free 
transitions, atomic free-free transitions, and molecular transitions. Currently, 
there are two codes available, which provide detailed information on these 
transitions. One code is NEQAIR developed by Park (Ref. 7) and another 
LORAN developed by Hartung (Ref. 8). Because NEQAIR code involves large 
database, they are rarely applied in practical problems. In this study, LORAN 
code is used to calculate the air plasma radiative properties. 

With the determination of k w and j' a , the RTE can be solved by either 

deterministic or stochastic approach. Unlike the most CFD equations, RTE is 
integral differential equation and numerical treatment is, thus, different from CFD 
approach. Currently, there are several methods available for solving non- 
equilibrium radiative heat transfer which include Monte Carlo method, P-1 
method, quadromoment method, discrete ordinates method (DOM), discrete 
transfer method (DTM), etc. The Monte Carlo is accurate but too time-consuming 
for practical applications. The P-1 method and quadromoment method are only 
accurate for optically thick medium. The DOM and DTM are mathematically very 
simple and they can provide accurate results for all optical ranges if the discrete 
direction number is reasonable large. Thus, the DOM and DTM are very suitable 
for modeling radiation with a participating medium. For the laser propelled vehicle 
application, the high temperature region where non-equilibrium radiation is 
prominent is usually relatively small compared to the flow region considered. If 
the DOM is applied, RTE must be solved in the entire flow region, thus, much of 
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the CPU time is wasted in the regions where radiation is not important. However, 
with the use of the DTM, RTE can be solved in the designated region. Therefore, 
the DTM is selected to solve the RTE in this study due to its high efficiency in 
solving the RTE. 

4.2 Results of Non-equilibrium Radiation Code Testing 

Based on the analysis described above, a computer code has been 

developed which is able to simulate multi-dimensional radiative heat transfer on 
unstructured grids in laser induced air plasmas. Before this code was coupled 
into the UNIC-UNS code, the LORAN code has been tested in order to gain a 
confidence on its usage. The LORAN code is available from NASA Langley 
Research Center and it consists of over fifty subroutines. The original code was 
designed for coupling with ID flow problems in a CRAY computer and its 
execution includes two separate steps involving spectrum setting up and 
radiation calculation. For the convenience of the present application, the code 
structure of the LORAN has been modified significantly. The modified code can 
be used to any UNIX computers and can be easily coupled into general CFD 
codes. In order to validate the accuracy of the modified LORAN code, two 
problems were selected and calculated absorption coefficient and/or emission 
coefficient were compared with available other results. 

Figure 4 shows the predicted absorption coefficient and emission coefficient 
distributions for equilibrium air with P=2 atm and T=1 0,000 K. Compared to the 
corresponding results predicted by the original LORAN code presented in Ref. 9, 
solutions from two codes are essentially identical. Figure 5 shows the predicted 
absorption coefficient distribution for equilibrium air with P=1 atm and T=1 4,000 
K. These results are very close to the other solutions presented on page 669 of 
Ref. 10. It is noted that the atomic line contributions are excluded in Figure 5 in 
order to create the same conditions as used in Ref. 10. 
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With the validation of the modified LORAN code, the radiation code 
developed in this study was then coupled into the UNIC-UNS code to investigate 
the effects of plasma radiation on the flow field. Figure 6 demonstrates the global 
laser absorption efficiency, temperature, pressure variations for the cases without 
and with consideration of plasma radiation during time Ips to 20 ps. The 
temperature and pressure are presented for a specific point closed to the focal 
point. Comparisons of results at different cases clearly indicate that plasma 
radiation has important effects on the flowfield at high temperatures, especially 
on the temperature distribution. Therefore, accurate prediction of plasma 
radiation becomes essential for high temperature regions. The drawback of 
coupling radiation into flow solver is its tremendous CPU time. Currently, the 
CPU time involving radiation coupling may be orders magnitude higher than that 
not involving radiation coupling. So, future development of non-equilibrium 
radiation model must focus on the improvement of computational efficiency. 



Figure 4. Absorption and emission coefficients for equilibrium air conditions. 
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Absorption Coefficient, 1/cm 



Frequency, (eV) 

Figure 5. Absorption coefficient at equilibrium air conditions 
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Figure 6. Effects of plasma radiation on global laser absorption coefficient, 

temperature, and pressure. 
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5. BENCHMARK TEST CASE STUDY 

To demonstrate the validity and effectiveness of the numerical model 
developed herein in predicting the performance data of a laser light craft, three 
test cases with different power levels were investigated in the present study. The 
test conditions of these cases are listed below. The laser ray tracing testing and 
time-sequence flowfield solutions using the grid system below are shown in the 
following pages. 

• Vehicle Configuration-A of Myrabo’s Flight Tests (Closed Air Inlet) 

• Pulsed Laser at 1 0 Hz and 30 ps Pulse Width 

• Laser Energy Level Up To 800 J (Test Points: 400 J, 600 J and 800 J) 

• Total Number of Mesh Elements = 26,318 

• Assuming Axisymmetric Model 

• Employs an automatic spark ignition scheme that the spark energy is 
switched off when the laser absorption efficiency reaches 15 percent. Also, 
only 15 percent of the beam power is used during the spark ignition. 

• Assuming 30 percent laser absorption where electron number density (related 
to the plasma frequency) reaches the resonance level, i.e. 2.51372E+23/m 3 . 

• With Ray T racing and Radiation Coupling 

• Testing the Predictions of Detonation Wave Propagation and Vehicle Thrust 
Calculation 



Figure 7. Mesh system used in the present laser light craft performance study. 
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5.1 Testing of the Ray Tracing Model 

The ray-tracing model is implemented in the UNIC-UNS code to simulate the 

laser ray reflection and focusing effects. At first, to avoid the inaccuracy caused 
by possible numerical noise in the wall curvature calculation, an analytical wall 
contour equation is used in the model. Later, a high accuracy method is 
employed to recover the analytical solution accurately based on the discrete grid 
points. The latter approach is general and produces accurate results. 


Figure 8 shows the results of ray tracing for the model-A configuration. A 
close up view near the focal point is shown in Figure 9. The red dots represent 
the ray-traced locations. The results give correct representation of the laser 
reflection and focusing effects. 



Figure 8. Model-A laser rays calculated by the ray-tracing model. 



Figure 9. Close up view near the focal point. 
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5.2 Flowfield Solutions 

The numerical computation is started with an assigned region of spark ignition 
near the laser focal point (a pie-shaped region covered by the beam). It is 
assumed that during this period, only 15 percent of the beam energy is converted 
to thermal energy. The spark ignition period is terminated when the laser 
absorption efficiency, calculated by the radiation model, researches 15 percent. 
At this time the maximum gas temperature reaches above 39,000 degrees K and 
the air plasma and a detonation wave are formed around the focal point. Local 
electron number density continues to grow until it reaches a level that the plasma 
frequency would be resonant with the laser frequency. When this level is 
reached, it is assumed that the laser energy is 100 percent absorbed by the 
electron and 30 percent of this energy is converted to thermal energy. 

Three cases with laser power level of 400 Joules, 600 Joules and 800 Joules 
are simulated. The vehicle thrust is calculated be integrating forces along the 
surfaces. Integration of the thrust with time is then recorded to calculate the 
coupling coefficient, which is defined as: 

J ( Vehicle - Thrust)dt 

Coupling Coefficient = ~z ; x 10 , in (dynes-sec/Joules) 

Other data, such as cumulative spark and laser energy and temperature, 
pressure and electron number density ranges are also recorded for analysis. 
Figures 10 to 14 shows the evolution of the detonation wave in pressure contours 
time sequence for the 400 Joules case. Temperature and electron number 
density contour plots at 19 microsecond are shown in Figure 15 to 17. Figure 15 
shows the case without the plasma frequency limit. In Figure 16, this limit is 
turned on, which shows that the laser beams are terminated at locations where 
the electron number density is on or above the resonance value. 
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Figure 10. Pressure distributions (ATM) at 5 microseconds. 



Figure 1 1 . Pressure distributions (ATM) at 10 microseconds. 
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Figure 12. Pressure distributions (ATM) at 15 microseconds. 



Figure 13. Pressure distributions (ATM) at 19 microseconds. 
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Figure 16. Temperature distributions (deg-K) at 19 microseconds with plasma 
frequency limit on the laser absorption. 



Figure 17. Electron number density distributions (1/m 3 ) at 19 microseconds. 
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5.3 Performance Data Comparisons 

The integrated coupling coefficient curves shown in Figure 18 reveal a typical 
initial rise due to the pressure waves acting on the vehicle surface. After 0.00017 
seconds, the coupling coefficients start to decrease due to the low pressure 
formed inside the cavity of the vehicle. Steady-state values of the coupling 
coefficients are obtained at around 0.001 1 second. For each case, the steady- 
state coupling coefficient is obtained in 2500 time steps. Figure 18 shows that the 
predicted steady-state coupling coefficients for the 400, 600 and 800 Joules 
power level cases are 13.8, 14.5 and 14.7 dyne-sec/Joules, respectively. These 
results are very close to the measured data of 12.6 dyne-sec/Joules. The data 
comparisons would have been better if the plasma radiation model is employed 
in the computation. However, it would take considerably large amount of CPU 
time to calculate just one case. Better scheme is needed to improve the 
computational efficiency of the plasma radiation model such that it can be 
included in the present analysis. 



Figure 18. Predicted coupling coefficient of the model A laser lightcraft. 
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6. CONCLUSIONS 


In the present project, a computational model has been developed for the 
prediction of laser propelled lightcraft performance. The present model involves 
the coupling of a time-accurate unstructured computational fluid dynamics (CFD) 
model and a robust radiative heat transfer model. Finite-rate chemistry model 
with high-temperature gas thermodynamics property is employed in the CFD 
code, UNIC-UNS. Electron and vibrational temperatures are also calculated in 
the present version of UNIC-UNS. Air chemistry system with 1 1 species and 51 
reactions is used in the present analysis. 

The present radiation model involves a ray-tracing model with a modified 
discrete transfer method for laser beam energy-absorption, propagation, 
reflection and focusing simulations. The modeling effort has also been spent on 
the investigation of the effects of high-temperature plasma radiation on the 
predicted plasma characteristics. It is found that the plasma radiation model does 
not have a major effect on the predicted strength of the detonation wave. 
However, it takes considerable amount of computer time to produce an answer. 
Therefore, the plasma radiation model is not included in the final performance 
data calculations. 

Data comparisons of the predicted coupling coefficients for different power 
level have revealed good correlation with the measured data. The success of the 
present model is attributed to the models incorporated that reflect the physical 
processes occurred during the initial phase of the laser pulse. This directly 
affects the strength of the detonation waves that is formed through converting the 
laser energy to thermal energy. The faster this process occurs the stronger 
pressure it produces and thus higher thrust of the vehicle. Another limiting factor 
is the resonance of the plasma frequency with the laser frequency. This controls 
the duration that the laser beam can be focused effectively for thrust generation. 
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